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We report results of a systematic analysis of matter-wave gap solitons (GSs) in three-dimensional 
self-repulsive Bose-Einstein condensates (BECs) loaded into a combination of a cigar-shaped trap 
and axial optical-lattice (OL) potential. Basic cases of the strong, intermediate, and weak radial 
(transverse) confinement are considered, as well as settings with shallow and deep OL potentials. 
Only in the case of the shallow lattice combined with tight radial confinement, which actually 
has little relevance to realistic experimental conditions, does the usual one- dimensional (ID) cubic 
Gross-Pitaevskii equation (GPE) furnish a sufficiently accurate description of GSs. However, the 
effective ID equation with the nonpolynomial nonlinearity, derived in Ref. [Phys. Rev. A 77, 
013617 (2008)], provides for quite an accurate approximation for the GSs in all cases, including 
the situation with weak transverse confinement, when the soliton's shape includes a considerable 
contribution from higher-order transverse modes, in addition to the usual ground-state wave function 
of the respective harmonic oscillator. Both fundamental GSs and their multipeak bound states are 
considered. The stability is analyzed by means of systematic simulations. It is concluded that 
almost all the fundamental GSs are stable, while their bound states may be stable if the underlying 
OL potential is deep enough. 
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I. INTRODUCTION 

Matter-wave gap solitons (GSs) are localized modes 
that can be created in elongated Bose-Einstein conden- 
sates (BECs) loaded into one-dimensional (ID) optical- 
lattice (OL) potentials, with the intrinsic nonlinearity 
induced by repulsive interactions between atoms. GSs 
have been the topic of a large number of original works. 
Results produced by these studies were summarized in 
several reviews [IHll • Within the framework of the mean- 
field approximation, the description of matter-wave pat- 
terns is based on the Gross-Pitaevskii equation (GPE) 
for the macroscopic wave function ip [6|: 

a) 

which has proven to be very successful in reproducing 
experimental results for the zero-temperature BEG (for 
instance, for multiple vortices, as shown in detail in Refs. 

In this equation, M is the atomic mass, the norm 
of the wave function ■0 is unity, N is the number of atoms, 
and g — Airh^a/M is the interaction strength with a being 
the s-wave scattering length. In this work, we consider 
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only repulsive interactions (i.e., a > 0). Furthermore, 
VI (rj^) — {l/2)MLLi\r\ is the radial-confinement poten- 
tial and Vz{z) is the axial potential, which may include 
the axial harmonic trap and the ID OL, with depth Vq 
and period d: 



V,Xz) = (l/2)Ma;^z^ -j- Vb sin^ {-Kz/d) 



(2) 



The energy scale in the underlying (no-interaction) lin- 
ear problem is set by the OL's recoil energy. En = 
;i27r2/(2Md2). 

Most commonly, theoretical studies of GSs in elon- 

S ted geometries are carried out in terms of the ID GPE 



ih- 



dt 



'2M'dz^ 



+ V,iz)4> + giuN\(f>\^(t>, (3) 



where giD = 2ahioj^. This is an effective evolution equa- 
tion for the axial dynamics — described by the ID wave 
function 4'{z,t) — which can be derived from the full 
three-dimensional (3D) GPE after averaging out the ra- 
dial degrees of freedom under the assumption that the 
radial confinement is so tight that the transverse dynam- 
ics are frozen to zero-point oscillations. These condi- 
tions, however, are not easy to realize using typical ex- 
perimental parameters [H , and when such conditions are 
not met the transverse excitations can no longer be ne- 
glected, making an essentially 3D analysis necessary. In 
this work, we aim to investigate different physically rele- 
vant regimes and capture 3D effects in the generation and 
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stability of matter-wave GSs in elongated BECs. This 
analysis should make it possible to determine the range 
of validity of the ID GPE for the description of the GSs 
under typical experimental conditions, as well as charac- 
teristic features of the fundamental and higher-order GSs 
in realistic situations. 

In principle, when the transverse excitations are rel- 
evant, Eq. ([3]) fails and one has to resort to the full 



or, alterna- 



3D GPE H]), as recently done in Ref. 
tively, use extended ID models that, within the frame- 
work of certain assumptions, take into account effects of 
higher-order radial modes on the axial dynamics of the 
condensate, leading to effective ID equations with the 
cubic-quintic [l5| or nonpolynomial [l6l . [l7j nonlinear- 
ities. In this work, we will consider both the full 3D 
GPE and the effective ID model with a nonpolynomial 
nonlinearity, which was derived, for the case of the self- 
repulsive nonlinearity, in Ref. ■ The latter one, which 
represents a simple generalization of the usual ID GPE, 
that reduces to Eq. ([3]) in the appropriate limit, has 
demonstrated an excellent quantitative agreement with 
experimental observations [18l - [2l| (chiefly, for delocalized 
dark solitons, which are natural patterns in the case of 
the self-repulsion; however, the comparison was not re- 
ported before for localized GS modes). We will demon- 
strate that, while the range of applicability of the ID 
GPE ^ is severely limited in realistic situations, the 
above-mentioned generalization gives a good description 
of stationary matter-wave GSs in virtually all cases of 
practical interest. 

The paper is organized as follows. The model is formu- 
lated in Section II, where we also recapitulate the deriva- 
tion of the effective nonpolynomial ID equation following 
the lines of Ref. [l^l, as this derivation is essential for 
the presentation of the results for the GSs. The main 
findings are collected in Sections III, where families of 
GS solutions are reported in several physically relevant 
regimes for strong, intermediate, and weak transverse 
confinement and shallow or deep OL potential. Both 
fundamental GSs and their multipeak bound states are 
considered. Only in the case of tight confinement com- 
bined with a shallow OL does the ordinary ID GPE pro- 
vide for a sufficiently accurate description of the GSs in 
comparison with results of full 3D computations. On the 
other hand, the effective nonpolynomial ID equation pro- 
vides for good accuracy in all cases; for the fundamental 
solitons and their bound states alike. The stability of 
the GSs is studied in Section IV by means of systematic 
simulations of the evolution of perturbed solitons. The 
conclusion is that the fundamental GSs are stable (ex- 
cept in a narrow region close to the upper edge of the 
band gaps), even in the case of strong deviation from 
the usual ID description. Multisoliton bound states are 
stable if the OL potential is deep enough. Conclusions 
following from results of this work, including applicability 
limits for the mean-field approximation and ID approxi- 
mations, are formulated in Section V. 



II. THE MODEL 

The model that was developed in Ref. [Ol for a BEG 
in the absence of OL potentials resorted to the adiabatic 
approximation, to neglect correlations between the trans- 
verse and axial motions. This approximation assumes 
that the axial density varies slowly enough to allow the 
transverse wave function to follow these slow variations. 
Because the OL imposes a new spatial scale, which may 
be more restrictive, it is necessary to find out if the adi- 
abatic approximation remains valid in the presence of 
the OL. To this end, we will now briefly recapitulate the 
derivation of the effective ID equation based on this ap- 
proximation. 

The starting point is the ansatz based on the factorized 
3D wave function 



V'(r,t) = ip{rj_;ni{z,t))(l){z,t), 



(4) 



with ni{z,t) being the local condensate density per unit 
length characterizing the axial configuration: 

ni{z,t) = N f (fr^\^{r^,z,t)\^ ^ N\<j){z,t)\\ (5) 



To derive Eq. ([5]), we have assumed that the transverse 
wave function ip is normalized to unity. Next, the sub- 
stitution of Eq. into the 3D equation ([T]) leads to 



2M 



2M dz^ 



h? dip d(j) 
M dz dz ' 



(6) 



Because of the very different time scales of the axial and 
radial motions, it is natural to assume that the slow axial 
dynamics may be accurately described by averaging Eq. 
([6]) over the fast (radial) degrees of freedom. Doing this, 
one obtains 

(7) 



Hiini)^ d'^rj_ip*[- 



2M 



,V1 



2M dz^ 



+ V±+gni\tp\ ]tp 



(8) 

Since ni (z, t) enters the last term of Eq. (|8]) merely as an 
external parameter, it is clear that, whenever the axial 
kinetic energy associated with the transverse wave func- 
tion may be neglected; namely. 



2M dz^ 



(9) 
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the radial dynamics decouple and /ij_(ni) can be deter- 
mined without the knowledge of the axial wave function. 
Actually, the right-hand side of Eq. ([9|) is the chemi- 
cal potential of an axially homogeneous condensate char- 
acterized by the density per unit length ni and wave 
function ip. For a sufficiently small value of the linear 
density {ani <^ 1), the chemical potential of the lowest- 
energy state of this homogeneous condensate can be read- 
ily obtained perturbatively. In this case, to the lowest or- 
der, (p{r±) is given by the Gaussian wave function of the 
ground state of the radial harmonic oscillator, with the 
corresponding chemical potential being huj±{l + 2ani). 
As the linear density increases, more radial modes be- 
come excited and, in general, the ground state of the 
corresponding homogeneous condensate involves many 
harmonic-oscillator modes. 

In previous works, it was shown using different ap- 
proaches that, also in this case, an analytical solution 
can be constructed In particular, by using a varia- 
tional approach based on the direct minimization of the 
chemical-potential functional, it was shown that, for any 
(dimensionless) linear density ani, a very accurate esti- 
mate for the chemical potential of the ground state is 
given by huJ±^/T+'ianl [1^. This expression can be 
easily extended to incorporate the case in which the con- 
densate contains an axisymmetric vortex (TtI (see also 
Ref. [23], where the intrinsic vorticity was included 
into the derivation of the ID nonpolynomial nonlinear 
Schrodinger equation in the case of self- attraction) . How- 
ever, in this work we restrict the consideration to zero 
vorticity. Using Eq. we see that a sufficient condi- 
tion for the second derivative in z appearing in Eq. (|8]) 
to be negligible is 

K^[(p] < huj±^/l + 4:ani. (10) 

Taking into account an estimate, Kz[f] ^ fi-^/(2MA^), 
where is the characteristic length scale in the axial 
direction, Eq. (|10|) can be rewritten as 

« hw^Vl + ^am. (11) 

When this condition holds, /i^(ni) coincides, to a good 
approximation, with the transverse local chemical poten- 
tial of the stationary radial GPE: 

(^-^^l + VAr±) + gnilipf^ if ^ ^M^, (12) 
and is given by 

IJ'±{ni) = haj±\/l + 4ani. (13) 

Finally, substituting Eqs. (US]) and (O into Eq. ^ 
and taking into account that is a real normalized wave 
function (which implies the vanishing of the integral in 
the last term), we arrive at the following effective ID 
equation to govern the slow axial dynamics of the con- 
densate: 



We note that the contribution from interatomic in- 
teractions enters the above equation through term 
hl^J_{^/T^h4xlNy)\^ ~ 1)0, which vanishes for iV 0. 
Thus, Eq. ([T^ incorporates an energy shift fiuj±, which is 
irrelevant for the dynamics but simplifies the form of the 
equation and makes the global chemical potential that 
follows from this effective ID equation exactly coinciding 
with the corresponding 3D result. 

The above derivation demonstrates that the validity of 
Eq. (HH) relies on two conditions: 

(i) The typical time scale of the axial motion must 
be much larger than the typical time scale of the radial 
motion wj^). 

(ii) The axial kinetic energy Kz [tp] associated with the 
transverse wave function must be negligible. 

The former requirement is necessary to allow the ra- 
dial wave function to adiabatically follow the axial dy- 
namics. While the specific temporal scale At depends 
on the particular initial conditions, typically, in the pres- 
ence of an OL, the fulfillment of this condition gets more 
difficult as the lattice period d decreases, or the linear 
density ani increases. In particular, a sufficient condi- 
tion for the validity of the adiabatic approximation is: 

S> aj_, aril ^ 1, with a± = ^ h/{Mijj±) being the ra- 
dial harmonic-oscillator length. Nevertheless, such a con- 
straint, which is hard to satisfy in realistic situations, is 
not a necessary condition. Actually, for stationary states 
At oo and condition (i) always holds. In the present 
work we are interested in this case, which is the most 
relevant one to seek for matter-wave GSs. 

A sufficient condition for the fulfillment of condition 
(ii) is given by the inequality (jlip . In the presence of an 
OL, the characteristic length scale A^ is typically on the 
order of the lattice period d, hence Eq. ([TT1) becomes 

Eb. < fta;±%/TT4arH, (15) 

where Eji is the corresponding recoil energy. We will 
demonstrate that this condition is well satisfied for sta- 
tionary GSs in most cases of interest. 

It is clear that, when the linear density is small enough 
(aril *C 1), Eq- (El): which exhibits the nonpolyno- 
mial nonlinearity, reduces to the usual ID GPE ([3]) with 
the cubic nonlinearity. This is the quasi-lD mean-field 
regime, which corresponds to condensates whose radial 
state, as given by a solution to Eq. (fT2|) . may be well 
approximated by the Gaussian ground state of the cor- 
responding harmonic oscillator [OjlH]- Thus, Eq. (ITi)) 
represents an extension of the ID GPE for condensates 
with larger linear densities or, equivalently, larger mean- 
field interaction energies. In such condensates, the radial 
ground state satisfying Eq. ((T^ is, in general, a linear 
combination of many harmonic-oscillator modes (a simi- 
lar situation takes place in the derivation of the effective 
ID equation for fermionic gases by means of the density- 
functional approach [25|). 

Inspection of Eqs. ([T]) or ([Ti)) demonstrates that the 
dynamical problem is fully controlled by four dimension- 
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less parameters, which may be defined as 

Er n_ uj^ Na 

huj± ' Er ' Ll)±_ ' aj_ 



(16) 



As said above, the recoil energy Er sets the energy scale 
of the underlying linear problem, while the nonlinear cou- 
pling constant Na/aj_ determines the order of magnitude 
of the mean-field interaction energy in units of the radial 
quantum hujj_ . Note that the dimensionless linear density 
ani is proportional to Na/a±. 

In this work, we assume the axial confinement to be so 
weak that the corresponding harmonic-oscillator length, 
az = y^h/{Mujz), is much larger than the lattice period 
d. Under these conditions, it is safe to neglect the mod- 
ulation induced in the condensate density by the axial 
harmonic trap and set lJz = 0, so that we are left with a 
uniform ID lattice potential acting along the axial direc- 
tion. 



III. MATTER- WAVE GAP SOLITONS IN 
DIFFERENT PHYSICAL REGIMES 

The stationary solutions of the 3D GPE are wave func- 
tions of the form tp{r,t) = ■!/'o(r) exp(— j/xt), with 
obeying the time-independent GPE: 

f^^o = + VAr±) + Vz{z) + gN jV-ol') ^o, 

(17) 

where /i is the chemical potential. When Eq. (ITi)) is 
applicable, one can instead generate the stationary axial 
wave function 4>i){z) by solving the effective ID equation 



2M dz^ 



+ V4z) + huj±^l + 4aN\(j>o\'^^ 00, 



(18) 

which, in the limit of ani <^ 1, reduces to the time- 
independent version of the usual ID GPE 

A^'/'o = i- — -^ + Vziz)+giuN\(j>o\^ + hLj^_ 

(19) 

Matter-wave GSs are characterized by a chemical po- 
tential lying in a band gap of the energy spectrum of the 
underlying linear system. To construct solutions for real- 
istic 3D matter-wave GSs in the effectively ID geometry, 
we looked for numerical solutions to Eqs. (IT71) and (|18p in 
different physically relevant regimes, and compared the 
obtained results with those produced by the usual ID 
GPE (fT9| to determine to what extent 3D contributions 
are relevant. In particular, the comparison allows us to 
estimate the accuracy and range of applicability of the 
effective ID equations and (ITOl) . 

Experimentally relevant values for the OL period d 
range from 0.4 to 1.6 /im which implies that, for the 
condensate of ^^Rb, ER/{2-Kh) ranges from 3.6 kHz to 
220 Hz. Taking into account that, typically, uj±/2ti < 1 



kHz, we conclude that Er/Hu!± > 1/4. On the other 
hand, for the condensate of ^^Na one has Efi/huj± > 1 
p6| : therefore, in this case the GS energy is sufficiently 
large to excite higher modes of the radial confinement, 
which makes 3D contributions always relevant. Since the 
*^Rb condensate is most relevant for the experimental re- 
alization of GSs in the quasi-lD setting [23, H^, in what 
follows below we use particular parameters of this con- 
densate to illustrate our results. Nevertheless, the results 
of the present work, which are always expressed in terms 
of relevant dimensionless parameters, are valid for any 
BEG. Actually, this is a direct consequence of the scaling 
properties of Eqs. ©, ©, and (fT4)) . 



A. Tight radial confinement: iSjj/fejx <C 1 

In this case, the quantum of radial excitations is much 
greater than the typical energy scale in the linear prob- 
lem. Note that, to realize this regime in a ^^Rb conden- 
sate, even in an OL of period d ~ 1.6 /im, one needs to 
use a harmonic trap with radial frequency lo ±/2t: > 2400 
Hz. 

Figure [Ija) shows the dimensionless chemical poten- 
tial of the noninter acting 3D condensate with ER/huj± ~ 
1/10, 



^ EE (/i - 1tIO^_)/Er, 



(20) 



as a function of the dimensionless lattice depth, s = 
Vq/Er. Since in this work we are interested in GSs with 
zero vorticity, this diagram, which represents the band- 
gap structure of the underlying linear problem, has been 
obtained from the zero-vorticity solutions of the linear 
version of Eq. p7|) . Regions I, II, and HI correspond 
to the lowest finite band gaps, which separate shaded 
bands, where linear solutions exist. An important point 
is that, within the region of interest in the parameter 
space (for Vq/Er < 25 and up to the third band gap), 
this 3D diagram is indistinguishable from that obtained 
using the linear version of the effective ID equation ([T5)) 
[which, obviously, coincides with the linear version of the 
stationary ID GPE (HH)]. In fact, Eq. ^ leads to a 
band-gap diagram that differs from Fig. [Ija) solely in 
the band marked by the arrow, which does not appear in 
the ID case. This extra band is, essentially, a replica of 
the lowest one, shifted up in energy by 2haj±/ER. Tak- 
ing into account that the energy spectrum of the radial 
harmonic oscillator is given by 



E = {2nr 



(21) 



where Ur = 0,1,2,... is the radial quantum number and 
m = 0,±l,±2,...is the axial angular-momentum quan- 
tum number, it is evident that the additional up-shifted 
band corresponds to the first-excited radial mode with 
m = 0. In the notation of Ref. fl3|, it corresponds to 
quantum numbers (n = 1, m = 0, rir — 1), with n being 
the band index of the ID axial problem. Thus, the ap- 
pearance of this band is a purely 3D effect that cannot 
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FIG. 1: (Color online) (a) Band-gap structure of a noninter- 
acting 3D BEC with Er/IuiJx^ — 1/10, as a function of the 
dimensionless lattice depth s = Vq/Er. The representative 
cases s = 2 and 20, indicated by vertical dashed lines, are 
considered in more detail in (b) and (c), respectively, which 
show the dimensionless chemical potential /I as a function of 
the quasimomentum q in the first Brillouin zone (in units of 
Tv/d). The right panels display the location of the *^Rb gap 
solitons considered in this work (points A-G). 



be accounted for by the above ID models. Since the en- 
ergy shift increases as Efi/haj± decreases, it is clear that, 
within the parameter region of interest. Fig. [TJa) is uni- 
versal in the sense that it is valid (for both ID and 3D 
systems) for any Eii/hu!± ^ 1. 

In this work, we restrict ourselves to the case of (1, 0, 0) 
GSs; that is, the solitons bifurcating from the lowest- 
energy band. Three-dimensional matter- wave GSs with a 
nontrivial radial structure have been studied in Ref. [3] . 
In this subsection, we consider OLs with En/hui^ = 1/10 
and depth Vq/Er = 2 or 20. These two representative 
cases correspond to the dashed vertical lines in Fig. (Tfa) 
and are shown in more detail in Figs. [IJb) and HJc), 
which display the corresponding band-gap structure as 
a function of the quasimomentum in the first Brillouin 
zone (in units of ir/d). Bold red lines in Figs. (Hb) and 
[IJc) have been obtained from the linear version of the 
effective ID equation ([T5)) while thin black lines corre- 
spond to results obtained by means of the linear version 
of the 3D equation pT|) that cannot be reproduced with 
the ID model. As mentioned above, the only appreciable 
difference between the ID and 3D results comes from the 



contribution of the first-excited radial mode [see the top 
part of Fig. (He)]. The case of Vq/Er = 2, displayed 
in Fig. [TJb), corresponds to the condensate in relatively 
shallow OL potentials. Under these circumstances, the 
linear energy spectrum does not differ too much from 
that corresponding to the translationally invariant case 
iVz = 0), and the band-gap structure exhibits wide en- 
ergy bands separated by relatively small gaps, which are 
tangible only in the lowest part of the spectrum. On the 
contrary, for Vq/Er = 20 [the tight-binding regime, see 
Fig. HJc)] the condensate is trapped in the deep periodic 
potential. In this case, the lowest part of the linear spec- 
trum is dominated by large gaps separating relatively 
narrow energy bands. The panels on the right side of 
Figs. [IJb) and [IJc) show the location, with respect to 
the corresponding linear band-gap diagram, of the GSs 
that will be considered below (points A-G). The horizon- 
tal axes in these panels indicate the number of atoms in 
each GS for the *^Rb condensates (with the s-wave scat- 
tering length a = 5.29 nm) in the trap with radial fre- 
quency uj±/2-K = 2400 Hz. For other parameter values, 
the GS family is described by the same plots, if consid- 
ered in terms of the above-mentioned nonlinear coupling 
constant, Na/a± [a = 5.29 nm and uj±/2tt = 2400 Hz 
correspond to a/a± — 0.024]. In other words, the num- 
ber of atoms in a GS created in the condensate with 
scattering length a' and transverse confinement radius 
a'j_, which is tantamount to the GS with particular val- 
ues of a, a± and N, is given by 



N' 



-N. 



(22) 



GS solutions have been obtained by numerically solv- 
ing the full 3D equation P7|) . as well as the effective 
ID equations ([T^ and ([TO]) . To this end, we have im- 
plemented a Newton continuation method, based on a 
Laguerre-Fourier spectral basis, that uses the chemical 
potential fi as a continuation parameter. To ensure the 
convergence, methods of this kind require initiation of the 
iterative procedure with a sufficiently good initial guess. 
This means that one needs a good estimate for both the 
axial and the radial parts of the wave function. Our com- 
putations used the following initial ansatz for the wave 
function: 



Ta I 



: exp 



2r2c 



M^). (23) 



where the z-dependent condensate's width, expressed in 
units of a_L , is 



T={l + 2aN\Mz)\' 



>l/4 



(24) 



and (t)o{z) — ^Jk^jl sech(fcoz) is the axial wave function. 
The number of particles, iV, and fco are free adjustable 
parameters. Note that the radial part of the wave func- 
tion is the same as the variational solution of Eq. 
P^ . which was found in Ref. [l7l |. 
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FIG. 2: (Color online) Atom density of the gap soliton corre- 
sponding to point B in Fig. [TJb), displayed as (a) an isosurface 
taken at 5% of the maximum density and (b) as a color map 
along a cutting plane containing the z axis, (c) Dimension- 
less axial density an\ obtained from the 3D wave function, 
as prescribed by Eq. ((5| (open circles) along with the cor- 
responding predictions from the nonpolynomial ID equation 
(fT8)) (solid red line) and the ordinary ID GPE im (dashed 
blue line). 



1. Shallow optical lattice: Vo/Er = 2 

Point A in Fig. [TJb) corresponds to a fundamental GS 
located near the bottom edge of the first band gap. It 
looks qualitatively similar to that shown below in Fig. 
[3l but with the peak linear density ani(O) ~ 0.04. In 
this case, the effective ID equations ([T8l) and (fT9|) yield 
results which are indistinguishable, on the scale of the fig- 
ures, from those obtained from the full 3D equation ([T7|. 
This is not surprising because, for these parameters, con- 
dition ([Tsl) is certainly satisfied, and inequality ani ^ 1, 
which guarantees the validity of the ID GPE ([T9| . also 
holds. The (dimensionless) chemical potential of this GS 
is Jl = 1.75, which implies /i — 1.175 huj± ~ hence 
the radial wave function of the condensate should not 
differ too much from the ground state of the correspond- 
ing harmonic oscillator. These fundamental GSs, how- 
ever, can only accommodate 9 particles (for the ®^Rb 
condensate) , which is too small to use the mean-field ap- 
proximation. Yet these solutions play an important role 
as building blocks of higher-order GSs. Point B in Fig. 



FIG. 3: (Color online) Same as Fig. [2] but for the gap soliton 
corresponding to point C in Fig. [TJb). 



[Ijb) corresponds to one of these higher-order solitons, 
which is built as a symmetric linear combination of three 
fundamental GSs, see Fig. [21 Its chemical potential and 
number of particles are Jl = 1.75 and N = 35. Note that, 
for the relatively shallow OL {Vq/Eji = 2), interference 
effects arising from the overlap between fundamental GSs 
sitting in adjacent lattice cells play an important role. 
For this reason, the contrast between the three central 
peaks and minima separating them is rather low in Fig. 
m and the total number of particles in the compound 
soliton does not coincide with the sum of the numbers of 
its fundamental constituents. Extending this procedure, 
one can readily build a sequence of compound GSs with 
an increasing number of peaks. 

Panels (a) and (b) in Fig. [5] display the 3D density of 
the three-peak GSs as, respectively, an isosurface (taken 
at 5% of the maximum density) and a color map in the 
cross-section plane drawn through the z axis. From these 
results, which have been obtained by solving the full 3D 
equation (I17p . we have also derived the axial linear den- 
sity ni{z), defined as per Eq. ([5]). This density is shown 
in Fig. HJc) by open circles, along with the correspond- 
ing results obtained from the effective ID equations ([T5)) 
and The OL potential is also displayed by the dot- 
ted line (in arbitrary units). It is seen that the effective 
ID equation ^TE\i (solid red lines) reproduces the 3D re- 
sults very accurately. The ID GPE ([TO)) (dashed blue 
lines) gives a slight discrepancy (~ 3%) against the 3D 
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results. This discrepancy, which would be invisible in the 
isolated fundamental GSs that build the compound, may 
also be explained by effects of the interference between 
the constituents. 

Point C in Fig. [Ijb) indicates the location of a fun- 
damental gap soliton in the {N, Jl) plane, which sits near 
the top edge of the first band gap. It contains 28 particles 
and has chemical potential Jl = 2.42, which corresponds 
to fi = 1.2A2huj±. This quantity is again very similar 
to hu!±, so that one expects the ID GPE to be a good 
approximation in this case. Figure [3] shows the 3D den- 
sity and the axial linear density ani{z) of this GS. In 
Fig. EJc) one can see that the ID GPE (HH) (the dashed 
blue line) reproduces the 3D results obtained from the 
full GPE p7|) (shown by open circles) within a 5% devi- 
ation. The effective ID equation (fT8|) (whose results are 
represented by the solid red line) is again more accurate 
and reproduces the 3D results with an error < 1%. As 
before, this fundamental GS can be used to build multi- 
peak compounds. 

As one moves upward in the bandgap, the 3D effects 
become stronger, which is a consequence of the fact that 
the corresponding number of particles and, thus, the non- 
linear interaction energy increase. The above results, 
however, indicate that, for the tight radial confinement 
{ER/hujj_ < 1) and shallow OL {Vq/Er < 2), the spe- 
cific 3D effects may eventually be neglected. In fact, un- 
der these conditions the maximum number of particles 
that can be accommodated in a fundamental GS in the 
first band gap is so small that inequality ani <C 1 always 
holds. This means that both the linear and the non- 
linear energies remain much smaller than the quantum 
of the radial excitation {Ej^, anihio± <C huj^), hence no 
higher transverse modes are significantly excited. There- 
fore, the situation considered in this subsection may be 
categorized as the quasi-lD mean-field regime, in which 
the ID GPE ()19|) accurately describes the matter-wave 
GSs, provided that condition ^ 1 holds (otherwise, 
the mean-field approximation will be invalid). 



2. Deep optical lattice: Vq/Er = 20 

In terms of the underlying OL, this is the case of the 
tight-binding regime, Vq ^ E^. The nonlinear tight- 
binding regime is realized when the potential depth Vb is 
much larger than both the linear and nonlinear (mean- 
field) energies (i.e., Vq >• Eji, ariihuj^). Under these 
circumstances, the condensate density is highly localized 
near potential minima, making the overlap between den- 
sities associated with different wells negligible. Point D 
in the first band gap of Fig. [Ijc) corresponds to a fun- 
damental GS with Jl = 6.79 and A^ = 18. Since its peak 
axial density is ani(O) ~ 0.2, and Vo/En = 20 implies 
Vq = 2 huj±, this fundamental GS belongs to the nonlin- 
ear tight-binding regime. Figure 2] shows a compound GS 
built of three such fundamental solitons. It corresponds 
to point E in Fig. [Tfc). Its chemical potential, Jl = 6.79, 



(a) 



(b) 




FIG. 4: (Color online) Same as Fig. [2] but for the gap soliton 
corresponding to point E in Fig. [IJc). 



is the same as that of the corresponding fundamental 
GS, and its number of particles, A^ = 55, now coincides 
with the sum of the number in its constituents (note that 
we always approximate A^ to the nearest integer). As 
can be seen from Fig. |4l this compound soliton exhibits 
three well-separated identical peaks (BEG droplets) , each 
one being practically indistinguishable from the above- 
mentioned fundamental GS. Figure HJc) shows that the 
results obtained from the effective ID equation (|18l) (the 
solid red line) agree very well with those produced by the 
fuU 3D GPE (H?]) (open circles). In particular, the ID 
equation yields the particle number N = 56, which is very 
close to A^ = 55, as given by the 3D solution. Since A^ 
represents a measurement of the norm of the wave func- 
tion, it is clear that the error in the number of particles 
refiects the error in the corresponding wave functions. 
The ID GPE ([19]), corresponding to the dashed blue line 
in Fig. HJc), gives A^ = 49, which implies an error ^ 10%. 

Point F in Fig. [T](c) corresponds to a fundamental 
GS located near the top edge of the first band gap. 
This soliton, which is very similar to that shown in 
Fig. [5l contains 71 particles and has chemical potential 
Jl — 11.5, which implies fi = 2.15 hLu±. At this value 
of fi there is a significant probability of the excitation 
of higher levels in the radial-confinement potential, 
which is corroborated by the fact that the peak axial 
density of this GS is ani(O) ~ 0.7. Because inequality 
ani <C 1 does not hold in this case, one cannot expect 
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FIG. 5: (Color online) Same as Fig. [2] but for the gap soliton 
corresponding to point G in Fig. \Tlc) . 



the ordinary ID GPE (HH) to be valid. In fact, it 
yields the number of particles TV = 53, which implies 
an error greater than 25%. Nevertheless, the effective 
nonpolynomial ID equation (jl8|) remains valid and quite 
accurate. It produces N = 73, which corresponds to a 
maximum error of 2.8%. This is not surprising, since 
condition (fT5|) . which guarantees the validity of the 
nonpolynomial equation, is satisfied in this case. 

The solution pertaining to point G in Fig. [TJc) is qual- 
itatively similar. It corresponds to the fundamental GS 
located near the top edge of the second band gap, with 
chemical potential Ji = 17.2 and N = 182. Figure [5] 
shows the 3D density and axial linear density ani{z) of 
this BEC droplet. In this case, ani(O) ~ 1.3, indicating 
a large contribution of excited radial modes. As a conse- 
quence, the ordinary ID GPE (|T9)) [the dashed blue line 
in Fig. [5{c)] fails to reproduce the 3D results (open cir- 
cles). It predicts N — 119, which means the error exceeds 
34%. As seen from Fig. Eljc), the effective ID equation 
(fist (the solid red line) remains accurate enough in this 
case, too. In particular, it yields N = 187, the respective 
error being 2.8%. 

The above results imply that, with the deep OL, the 
number of particles that can be accommodated in a fun- 
damental GS can be large enough to make the 3D char- 
acter of the wave function essential. In fact, even for the 
tight radial confinement {En/huj^ <C 1), which is the 
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FIG. 6: (Color online) (a) Band-gap structure of a noninter- 
acting 3D BEC with Eii/hw±^ = 1/4, as a function of the 
dimensionless lattice depth a = Vq/Er. The representative 
cases s = 2 and 20, indicated by vertical dashed lines, are 
considered in more detail in (b) and (c), respectively, which 
show the dimensionless chemical potential /I as a function of 
the quasimomentum q in the first Brillouin zone (in units of 
Tv/d). The right panels display the location of the ^^Rb gap 
solitons considered in this work (points H-M). 

most favorable case for the applicability of the ID ap- 
proximation, the usual ID GPE with the deep OL 
potential cannot be reliably used beyond the first third 
of the first band gap. Since the validity of the mean 
field treatment requires N ^ 1, the usual GPE cannot 
be used close to the bottom edge of the band gap, ei- 
ther. It is thus clear that the range of validity of this 
standard ID equation is limited. This conclusion be- 
comes even more severe if one takes into account that 
the tight-radial-confinement regime is not easy to real- 
ize using typical experimental parameters. Our simula- 
tions indicate, however, that the effective nonpolynomial 
ID equation (fT5|) properly describes such essentially 3D 
situations, providing for an accurate description of the 
fundamental GSs in the entire band gaps. 

B. Intermediate radial confinement: Er/1ujJ±^ ~ 1/4 

This regime is of particular interest because it corre- 
sponds to typical experimental parameters. It can be 
realized, for instance, with a ®^Rb condensate in an OL 
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of period d = 1.55 /im, radially confined by a harmonic 
potential with lu±/2tt = 960 Hz. The particle numbers 
N of the GSs considered below are given for these physi- 
cal parameters, and, they can be converted into values of 
N for other situations by means of Eq. p2|) (now, with 
a/a± = 0.015) . 

Figure [6][a) shows the band-gap structure produced 
by the linearized version of Eq. (|17p for zero-vorticity 
modes. Recall that the linearization of the effective ID 
equations and pJl) yields, instead, the diagram dis- 
played in Fig. [TJa) (except for the narrow band marked 
by the arrow, which, as already mentioned, does not ap- 
pear in the ID approximation). It is seen that the ID 
approximation cannot reproduce the 3D picture resulting 
from the excitation of the radial modes, which manifest 
themselves in Fig. |6l[a) as replicas of bands shifted up 
by integer multiples of 2huj^/ Eji. Since this quantity is 
smaller in the present case than it was before, the ef- 
fect of these contributions in the region of interest is now 
stronger. While it is obvious that the ID equations can- 
not account for the 3D effects originating from the shifted 
bands, we will demonstrate that these equations may still 
produce an accurate description of common GSs (i.e., 
those originating from the m = Ur = energy bands). 

Figures [SKb) andHJc) display the band-gap structure 
as a function of the quasimomentum for lattice depths 
Vq/Eh = 2 and 20, respectively. Thin black lines in these 
figures represent 3D results that cannot be reproduced 
by the ID models, and points H-M in (N, ]!) plane mark 
examples of GSs that will be considered below. 



1. Shallow optical lattice: Vo/Er = 2 

Point H in Fig. [IKb) corresponds to a fundamental GS 
with /i = 2.19 and iV = 50, which is similar to the one in 
Fig. [21 but with the peak axial density ani{Q) ~ 0.2. In 
this case, the effective ID equation (|18l) predicts = 51, 
thus corresponding to a 2% error, while the ID GPE 
predicts = 45 (an 11% error). Symmetric combi- 
nations of five such fundamental GSs generate the com- 
pound GS shown in Fig. [7l which contains 258 particles 
and corresponds to point I in Fig. ^h). The 3D den- 
sity of this soliton exhibits five weakly separated peaks, 
as shown in Figs. Efa) and[7fb) by means of the isosur- 
face and the color map in the cross-section plane drawn 
through the z axis. As seen in Fig. EKc), in this case 
the effective ID equation yields an error of 1.5% 
in the norm of the wave function {N = 262) while the 
use the ordinary ID GPE p9| generates an error of 12% 
{N — 226), showing that, already for this GS, the 3D 
effects play an important role. 

Point J in Fig. HKb) corresponds to a fundamental GS 
containing 75 particles and located near the top edge of 
the first band gap. It is qualitatively similar to the one 
shown in Fig. [31 but having /I = 2.42 and ani(O) ~ 0.26. 
In this case, Eqs. (fT5)) and yield results with an error 
~ 1% {N = 76) and 12% {N = 66), respectively 
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FIG. 7: (Color online) Atom density of the gap soliton corre- 
sponding to point I in Fig. [Hlb), displayed as (a) an isosurface 
taken at 5% of the maximum density and (b) as a color map 
along a cutting plane containing the z axis, (c) Dimension- 
less axial density an\ obtained from the 3D wave function, 
as prescribed by Eq. ((5} (open circles) along with the cor- 
responding predictions from the nonpolynomial ID equation 
p8)) (solid red line) and the ordinary ID GPE IHl) (dashed 
blue line). 



2. Deep optical lattice: Vq/Er = 20 

Points K, L, and M in Fig. ^c) correspond to funda- 
mental GSs in the case of the tight-binding underlying 
linear structure. The first soliton, which contains 19 par- 
ticles, with Jl = 5.33 and axial density ani(O) ~ 0.23, also 
corresponds to the nonlinear tight-binding regime and is 
thus highly localized at a single lattice site. This follows 
from the fact that, for Vo/Efi^ — 20 and Eft/huj± = 1/4, 
one has Vq/ huj± — 5, which is much greater than the 
above-mentioned value of ani(O). For this GS, located 
near the bottom edge of the first band gap, the ordinary 
ID GPE p^ yields an error > 11%, which continues to 
grow as one moves upward in the band gap. However, the 
effective nonpolynomial ID equation (|18p remains accu- 
rate within a 2.5% deviation. 

Point L in Fig. [SJc) indicates the position of a funda- 
mental GS near the top edge of the first band gap, with 
Jl = 11.5, N — 243, and ani(O) ~ 2.3. Since the inequal- 
ity aril <^ 1 does not hold in this case, the ordinary ID 
GPE (UnD is invahd, giving a 45% error {N = 133). Note 
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FIG. 8: (Color online) Same as Fig. [7]but for the gap soliton 
corresponding to point M in Fig. [SJ^c) . 



that, even though point L is relatively close to a 3D en- 
ergy band [the thin black line in Fig. El^c)], which cannot 
be reproduced by the effective ID equation p^ . it can, 
however, reproduce this GS within a 4% error {N = 253). 
The same is true for the fundamental GS displayed in 
Fig. m which corresponds to point M in Fig. IHc). It 
represents a BEG droplet containing 696 particles, with 
Jl = 17.2 and the peak axial density ani(O) ~ 5.4. Since 
we have Vq/ Huj± < ani(O) in the present case, the sys- 
tem is no longer in the nonlinear tight-binding regime. 
This is seen in Fig. [SJ where the density is not strongly 
localized around the minimum of the potential well. As 
seen in Fig. [H^c), the nonpolynomial ID equation p8)) 
yields results (the solid red line) that agree with the 3D 
picture produced by the full GPE ([T7]) (open circles), 
within a 3% deviation [Eq. ([TH]) yields N = 719 in this 
case]. On the contrary, the ID GPE (|19l) . which gives 
results (the dashed blue line) with an error > 55% (it 
yields N = 297), clearly is not valid in this case. 

These findings demonstrate that, in the physically 
relevant regime of the intermediate radial confinement 
{Eji/tiujj_ > 1/4), even for the shallow OL the 3D effects 
may be important, and thus the usual ID GPE ([T5|) fails 
to reproduce correctly the axial density ani{z) and the 
particle content N. For Eji/huj± = 1/4 and potential 
depth Vq/Eh = 2, the fundamental GSs located in the 
first band gap, as predicted by the ID GPE equation, fea- 
ture the error ~ 12%, which is still larger for compound 



FIG. 9: (Color online) (a) Band-gap structure of a nonin- 
teracting 3D BEG with Eii/tiu± — 1, as a function of the 
dimensionless lattice depth s = Vo/Er- The representative 
cases s — 2 and 20, indicated by vertical dashed lines, are 
considered in more detail in (b) and (c), respectively, which 
show the dimensionless chemical potential /I as a function of 
the quasimomentum q in the first Brillouin zone (in units of 
Tv/d). The right panels display the location of the *^Rb gap 
solitons considered in this work (points P-R). 



GSs. For the deep OL {Vq/Er = 20), the ID GPE is 
valid only in a small region near the bottom edge of the 
first band gap. In general, this equation is applicable only 
where both conditions ani <^ 1 and N ^ 1 hold simul- 
taneously. As Eji/huj± increases, the relative strength of 
the radial confinement decreases, and the range of valid- 
ity of the ID GPE becomes more and more narrow. From 
the experimental viewpoint, the increase of Eii/huj± may 
be relevant because, in this way, one can easily increase 
the number of particles in the fundamental GSs. How- 
ever, the increase of EBjfvjj± also implies a decrease in 
the relative size of the radial-excitation energy quantum, 
which can compromise the stability of the solitons be- 
cause they can decay by exciting higher radial modes, 
even for condensates with a relatively small number of 
particles. We briefly consider this case below. The sta- 
bility properties of the GSs will be analyzed in Sec. IV. 
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C. Weak radial confinement: E[i/hw±_ > 1 

In this regime, the typical energy scale in the under- 
lying hnear problem is sufficiently large to easily excite 
higher transverse modes. As a representative example, 
we consider the case of Efj/hLu± = 1, which can be re- 
alized in a ^^Rb condensate in an OL with d = 1.55 
^m and radial-confinement strength uj±/2tt — 240 Hz. 
The particle contents of the GSs considered below corre- 
spond to such condensates. As before, the correspond- 
ing values of N can be converted into those correspond- 
ing to other situations by means of Eq. (this time, 
a/ajL = 7.6 x IQ-^). 

Figure [SJa) displays the band-gap structure of the un- 
derlying 3D linear problem, to be compared with Fig. 
[ija) which shows the band-gap structure obtained from 
the linearization of ID equations or ([T^. As before, 
the latter equations cannot account for the 3D contribu- 
tions from the excited radial modes, which account for 
the series of shifted bands separated by gaps of value 
2huj±/ER in Fig. M^a). Figures H^b) andlijc) display the 
linear band-gap structure as a function of the quasimo- 
mentum for Vo/Er = 2 and 20, respectively. 

Point P in Fig. [9jb) corresponds to a fundamental GS 
with Ji = 1.59 and N = 74. Its 3D density plot, shown 
in Fig. [TUl demonstrates that this fundamental soliton is 
spread over several lattice sites, which is a consequence 
of the fact that its chemical potential is very close to the 
first linear energy band, where only extended solutions 
of the stationary GPE exist. Figure fTUrc) shows that 
the results obtained from the effective ID equation 
(the solid red line) coincide with the 3D results (open 
circles) within 1% (it predicts N = 75). However, the 
ID GPE (the dashed blue line) gives rise to an error 
~ 10% (it predicts N = 66). We thus conclude that, in 
the weak-radial-confinement regime {Eii/huj± > 1), the 
latter equation is only applicable in a narrow region close 
to the bottom edge of the first band gap, even for shal- 
low OLs. In this regime, the 3D contributions play an 
important role in most cases. These contributions are, 
however, well accounted for by the effective ID equation 
(IT51) . For instance, for point Q in Fig. MJ^) [which corre- 
sponds to a fundamental GS near the top edge of the first 
band gap, with Jl = 2.42, N = 400, and ani(O) ~ 1.5] 
the ID GPE (dH) gives an error of 34% (N ^ 265), while 
the nonpolynomial ID equation (|18p limits the error to 
3.5% (AT = 414). Point R in Fig. ^c) is an example of a 
fundamental GS in a deep OL. It corresponds to a disk- 
shaped BEG droplet trapped in a single lattice cell [see 
Fig. [TlTb')]. which contains 2323 *^Rb atoms, and has 
/i = 11.5. Its axial linear density is qualitatively similar 
to that shown in Fig. [SJc), but with a maximum value 
arii(O) ~ 23. In this case, the ID GPE predicts the 
particle content N = 534, while the effective ID equa- 
tion (|18p yields N — 2437, which corresponds to an error 
~ 5% in the norm of the wave function. Thus, the ID 
nonpolynomial equation provides for a sufficiently good 
description of the stationary GSs even in the highly non- 
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FIG. 10: (Color online) Atom density of the gap soliton corre- 
sponding to point P in Fig. H^b), displayed as (a) an isosurface 
taken at 5% of the maximum density and (b) as a color map 
along a cutting plane containing the z axis, (c) Dimension- 
less axial density ani obtained from the 3D wave function, 
as prescribed by Eq. ((5} (open circles) along with the cor- 
responding predictions from the nonpolynomial ID equation 
(O (solid red line) and the ordinary ID GPE (dashed 
blue line). 



linear (ani 3> 1) weak-radial-confinement regime. 



IV. STABILITY ANALYSIS 

We have investigated the stability of the GS solutions 
by monitoring their long-time behavior after the applica- 
tion of a random perturbation [l^, [1^ . Specifically, we 
have perturbed the corresponding wave functions with a 
small-amplitude (^ 2%) additive Gaussian white noise 
and have monitored the subsequent nonlinear evolution 
for 1 s. To this end we have numerically solved the full 3D 
GPE ([T]), as well as the effective ID equation (|14p . using 
a Laguerre-Fourier pseudospectral method with a third- 
order Adams-Bashforth time-marching scheme. Numeri- 
cal integration of the 3D GPE for such a long time is a 
very demanding computational task. While it is possible 
for a certain GS to be metastable and decay on a time 
scale still longer than 1 s, in practice this time is long 
enough in comparison with the lifetime of the conden- 
sate per se, hence the soliton surviving for 1 s may be 
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FIG. 11: (Color online) Evolution in time, after the applica- 
tion of a 2% white-noise perturbation, of the fundamental gap 
solitons corresponding (a) to point P in Fig. [9jb) and (b) to 
point R in Fig. IHc). 



FIG. 12: (Color online) (a) Evolution in time, after the ap- 
plication of a 2% white-noise perturbation, of the compound 
gap soliton corresponding to point B in Fig. [ijb). Panel (b) 
displays the long-time behavior, after the application of the 
perturbation, of a stable three-peak soliton (see text). 



categorized as stable. 

We have found that both the full 3D GPE ^ and 
the effective nonpolynomial ID equation (ITil) lead to the 
same conclusions regarding the stability of the GSs. How- 
ever, once an unstable soliton begins to decay, one cannot 
expect, in general, the latter equation to reproduce the 
dynamical behavior correctly. The reason for this prob- 
lem is that, as already mentioned in Sec. II, in the pres- 
ence of the OL potential the above-mentioned adiabatic 
condition (i) is hard to fulfill in time-dependent settings. 
Of course, the same is true for the ID equation ([3]), which 
is a limit form of Eq. (jl4p and, consequently, has a much 
narrower range of applicability. The results presented be- 
low have been obtained using the full 3D GPE ([Ij. For 
each GS we have used at least two different basis sets 
and time steps to check that the results do not depend 
on details of the numerical procedure. 

Our simulations demonstrate that, in general, (1,0,0) 
fundamental GSs are stable, except in a narrow region 
close to the top edge of the band gaps. In particular, the 
solitons corresponding to points A, D, K, and P in Figs. 
[IJ \6\ and |9] remain stable up to t = 1 s. On the contrary, 
GSs in a shallow OL {Vo/Eji = 2), located close to the 
top edge of a band gap, such as those corresponding to 
points C, H, J, and Q in Figs. [U [U and 13 turn out to 
be unstable. This instability manifests itself as a steady 
decay of the norm of the GS through emission of radia- 



tion, on a time scale that increases as one moves away 
from the top edge of the band gap. In this regard, one 
finds that the GS corresponding to point H decays much 
slower than the other ones. As the lattice depth Vq/Eh 
increases, in general, GSs become more stable and the 
instability region shrinks. In particular, our simulations 
indicate that GSs such as those corresponding to points 
F, G, L, M and R in Figs. [Hil and |9] (which are funda- 
mental GSs in the deep OL, with Vq/Eh = 20, located 
close to the top edge of the band gap) also remain stable 
up to t = 1 s. This is in good agreement with previous 
analyses carried out in the context of deep optical lattices 
in terms of the ordinary ID GPE © (i,|l3. 

An important result is that (1,0,0) fundamental GSs 
remain stable even in the weak-radial-confinement regime 
{Eii/Huj± ~ 1). In this regime, GSs always have sufficient 
energy to excite higher radial modes, and, as a conse- 
quence, the 3D effects are always relevant and the usual 
ID GPE ([3]) fails. Figure [rTTa) displays the long-time be- 
havior, after the application of the perturbation, of the 
fundamental GS shown in Fig. [10] [it corresponds to point 
P in Fig. IHIb)]. This figure shows (by means of a color 
map) the evolution of the axial density ani{z,t), which 
has been obtained from the 3D wave function 4^(r±, z, t) 
by integrating out the radial dependence, as per Eq. ([5]) . 
The left panel in this figure shows the 3D condensate den- 
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sity at t = (i.e., just after the application of the pertur- 
bation) as an isosurface taken at 5% of the maximum den- 
sity, while the right panel represents the density at < = 1 
s. This GS has Ji — 1.59, which implies /i = 2.59huj±. 
Note that, despite the fact that this chemical potential 
is greater than the quantum hio± , the GS remains stable 
up to 1 s. The same is true for the fundamental GS cor- 
responding to point R in Fig. IHtc), which corresponds 
to a disk-shaped BEC containing more than 2300 *^Rb 
atoms with chemical potential fj, = 12.5huj± ^ hujj_. As 
Fig. [TlTb) shows, this GS also remains stable. 

Our simulations also demonstrate that the compound 
GS from Fig. [2l corresponding to point B in Fig. [TJb), 
is unstable. This is seen in Fig. Il^ a'). which shows how 
the soliton decays on a time scale ~ 120 ms after the 
addition of a 2% white-noise perturbation. This com- 
pound GS, which has Jl = 1.75 and = 35, was built in 
the shallow OL {Vq/Ej^ — 2) as the symmetric combina- 
tion of three fundamental constituents. We have found 
that, in such shallow lattices, GSs of this type are al- 
ways unstable. However, it is not difhcult to find stable 
solitons of this type in somewhat deeper lattices. An ex- 
ample is shown in Fig. I12f b'). which represents a stable 
three-peak soliton with Jl — 3.13 and A^ = 59, trapped 
in the OL with Vq/Er = 4 and EB./huj± = 1/10. Sim- 
ilar dynamics is observed for the GS displayed in Fig. 
[3 which corresponds to point I in Fig. IH^b). This 
is a five-peak compound generated in the shallow OL 
{Vq/Eh = 2) in the regime of the intermediate radial- 
confinement strength {Ejj/huj± ~ 1/4). Our simulations 
indicate that this compound soliton is unstable. In fact, 
no stable five-peak solitons were found in such a shallow 
lattice. On the other hand, it is not difficult to find stable 
compounds of this kind for deeper OLs. An example is 
the five-peak pattern with Vq/Eh = 4, Eii/huj± = 1/4, 
Ji — 2.83, and A^ = 212. In general, GSs naturally be- 
come more stable as the lattice depth increases. For 
deep OLs {Vo/Efi = 20), the stability of compound GSs 
is essentially identical to that of their fundamental con- 
stituents. The GS in Fig. |4l which corresponds to point 
E in Fig. [Tfc) , is an example of a stable three-peak soli- 
ton realized in a deep OL. 

V. CONCLUSIONS 

To appraise the physical relevance of the results re- 
ported in this work, it is pertinent to recall that the 
mean- field treatment of the GSs (gap solitons), based on 
the GPE, is valid if A^ ^ 1 and the condensate is suffi- 



ciently dilute so that a^n ^ 1, where a is the scattering 
length of the inter-atomic collisions and n is the atomic 
density. In shallow OLs (optical lattices) such conditions 
can be easily met, though the former one imposes a seri- 
ous limitation on the range of applicability of the usual 
ID GPE (fT9|) . which fails as N increases. In deep OLs, 
where the tunneling between adjacent cells is strongly 
suppressed, the above conditions must be fulfilled at each 
site of the lattice. Under these circumstances, GSs may 
have a rather high local density. A good estimate for the 
3D peak density n(0) = N\ip{0)\'^ in terms of the peak 
axial density, ni(0), can be obtained from the following 
relation [13] : 



■Ka]_^y^ + 4ani(0) 

Substituting the values of ani(O) obtained in Section III 
for different GSs into Eq. (l25t . one can easily verify 
that condition a^n{0) <C 1 is satisfied in all cases, which 
justifies the use of the mean-field description. The most 
critical situation occurs for the GS corresponding to point 
G in Fig. IHc), which has a^n{0) ~ lO"'^. For the GS 
corresponding to point R in Fig. |9l[c), one finds a^n{0) ~ 
5 X 10-5. 

Despite the fact that applicability conditions for the 
mean-field treatment can be readily met, it is much 
harder to justify the validity of the usual ID GPE ([3|), 
which requires both conditions arii <^ 1 and N ^ 1 
to hold simultaneously. Only in this case may the 3D 
effects be safety neglected. In most experimentally rele- 
vant situations, these conditions are not met, hence the 
applicability range of Eq. ^ turns out to be very lim- 
ited. On the contrary, it has been shown in this work 
that the effective ID GPE (fTi|) with the nonpolynomial 
nonlinearity provides for an accurate description of the 
stationary matter-wave GSs in most cases of practical 
interest. 

A relevant extension of the present analysis may be 
to GSs with intrinsic vorticity, as well as to mobility of 
stable solitons. It is also interesting to consider two- 
component GSs in the realistic 3D setting. 
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